###### Appendix: Table A4
###### Evidence of utility gains from punishment
gc(); rm(list = ls()); set.seed(12345)
setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) # Note: if you are not using R Studio this command will not work, set WD to source file location manually

############### SCRIPT SUMMARY #####################

packages <- c("dplyr", "ggplot2", "texreg")
lapply(packages, require, character.only = T)
source("functions.R")

### Load Datasets

dfMerge <- read.csv("data/cleaned/uganda_jun17_cleaned.csv", header=T, stringsAsFactors =F)
dfMerge <- dfMerge[!dfMerge$game_id_exp %in% c("Basic Aid", "Aid Ownership"), ]

########################################## MISC PREP AND VARGEN ######################################

## Re-level treatment variable according to which one you want
print(table(dfMerge$treatment_pooled))
newLabs <- c("Hidden VAT", "Direct Tax", "Visible VAT", "Windfall") 
origLabs <- c("Strong Hidden VAT", "Direct Tax", "Weak Hidden VAT", "Windfall")
treatCor <- setNames(newLabs, origLabs)
dfMerge$treatment <- treatCor[match(dfMerge$treatment_pooled, 
                                    names(treatCor))] 
table(dfMerge$treatment)
dfExp <- dfMerge
dfExp <- dfExp[dfExp$role == "C" & dfExp$round != 1 & 
                 !is.na(dfExp$treatment), ] # First round is a practice round

## Determine which columns to take from the model frames

colsInc <- c("outcome", "comparison", "estimate", "std.error", 
             "conf.low", "conf.high", "p.value", "n")

### Do expressive benefits from punishment exist?

require(estimatr)

refCat <- c("Visible VAT", 
            "Hidden VAT", 
            "Direct Tax", 
            "Windfall")

tabVars <- c("as.factor(punish)1", 
             "treatmentDirect Tax", 
             "treatmentHidden VAT", 
             "treatmentVisible VAT", 
             "prev_transfer", 
             "sbj_thresh", 
             "leader_transfer", 
             "loss_degree")


for(i in 1:length(refCat)){
  
  depVar <- c("l4_l3")
  dfUsing <- dfExp 
  dfUsing$treatment <- relevel(factor(dfUsing$treatment), 
                               ref = refCat[i])
  expVars <- c("as.factor(punish)", "treatment", "as.factor(enum_id)", 
               "item", "as.factor(round)", "prev_transfer", "sbj_thresh", "leader_transfer")
  controlVars <- c("age", "edu", "poverty_avg", "services_avg", "female")
  
  modEq <- as.formula(make_equation(depVar, c(expVars, controlVars)))
  modFit <- lm_robust(modEq, data = dfUsing, 
                      clusters = sbj_unique,
                      se_type = "CR2")
  
  dfMod <- tidy(modFit)
  dfMod <- dfMod %>% 
    filter(term %in% c(refCat, tabVars))
  dfMod$n <- modFit$nobs
}

### Does the degree to which you experience loss affect your benefit from punishment?

for(i in 1:length(refCat)){
  
  depVar <- c("l4_l3")
  dfUsingloss <- dfExp %>% 
    filter(punish == 1)
  dfUsingloss$loss_degree <- 10 - dfUsingloss$l2
  dfUsingloss$treatment <- relevel(factor(dfUsingloss$treatment), 
                               ref = refCat[i])
  expVarsloss <- c("loss_degree", "treatment", "as.factor(enum_id)", 
               "item", "as.factor(round)", "prev_transfer", "sbj_thresh", "leader_transfer")
  controlVars <- c("age", "edu", "poverty_avg", "services_avg", "female")
  
  modEqloss <- as.formula(make_equation(depVar, c(expVarsloss, controlVars)))
  modFitloss <- lm_robust(modEqloss, data = dfUsingloss, 
                      clusters = sbj_unique,
                      se_type = "CR2")
  
  dfModloss <- tidy(modFitloss)
  dfModloss <- dfModloss %>% 
    filter(term %in% c(refCat, tabVars))
  dfModloss$n <- modFitloss$nobs
}

dfAll <- as.data.frame(rbind(dfMod, dfModloss))
roundVars <- c("estimate", "std.error", "p.value")
dfAll[, roundVars] <- round(dfAll[, roundVars], digits = 2)

useCols <- c("term", "estimate", 
             "std.error", "p.value", "n")

(dfComb <- dfAll[, useCols])



